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Applications  of  the  Finite-Element  Method  to  the 
Problem  of  Heat  Transfer  in  a  Freezing  Shaft  Wall 


FU  LIANDI 


INTRODUCTION 

Artificial  freezing  has  been  widely  used  in  the  construction  of  trenches,  tunnels,  bridge  open¬ 
ings  and  other  underground  structures.  The  essential  problems  with  this  technique  are  the  deter¬ 
mination  of  temperature  change  and  the  control  of  the  freezing-thawing  process  in  the  develop¬ 
ing  frozen  wall.  The  use  of  this  technique  is  closely  related  to  the  engineering  design,  construc¬ 
tion  procedures,  economic  benefits  and  work  safety. 

In  the  area  of  Liang-Huai  Coal  Mine  in  China,  the  depth  of  the  shaft  is  generally  300  to  700  m. 
For  each  shaft,  about  40  freeze  pipes  are  placed  symmetrically  along  the  circumference  of  radius 
R.  The  depth  of  the  freeze  pipe  is  usually  identical  to  that  of  the  shaft.  The  freeze  pipes  may 
also  be  staggered  along  two  concentric  circles,  called  freeze-pipe-circles,  with  radii  /?  ]  and  /?2 . 

The  coolant  used  for  freezing  is  a  brine  solution  with  a  temperature  of  approximately  -30°C. 

By  circulating  the  cold  brine  in  pipes  and  the  subsequent  heat  extraction  between  the  pipe  walls 
and  soil,  a  frozen  wall  is  thus  formed.  The  process  of  frozen  wall  formation  may  be  divided  into 
four  stages  (Beijing  Institute  1975). 

In  the  first  stage,  which  is  called  precircumscribed  circling,  freezing  progresses  radially  at  each 
freeze  pipe,  and  a  thin  frozen  wall  will  appear  around  each  freeze  pipe.  After  a  certain  duration 
that  depends  on  the  distance  between  every  two  adjacent  freeze  pipes  and  the  specific  flow  con¬ 
dition  of  cold  brine,  all  the  small  frozen  walls  will  be  connected  to  each  other,  and  a  global, 
nearly  circular  frozen  wall  will  be  formed. 

After  the  formation  of  the  global  wall  in  the  first  stage,  the  frozen  wall  will  concurrently  grow 
inward  and  outward  until  its  thickness  and  temperature  reach  desired  values.  This  is  called  the 
active  freezing  stage. 

During  the  third  stage,  continuous  freezing  is  necessary  to  maintain  an  optimum  thickness  and 
temperature  of  the  frozen  wall  while  excavation  and  cement  lining  of  the  mine  shaft  are  taking 
place.  In  the  last  stage,  the  circulation  of  cold  brine  is  stopped  as  soon  as  the  shaft  work  is  com¬ 
pleted.  This  is  called  the  natural  recovery  period.  Of  the  four  stages,  the  active  freezing  stage  is 
most  important  for  predesigning.  This  paper  will  focus  attention  on  this  stage  and  try  to  calculate 
the  change  in  the  size  of  the  frozen  wall  with  the  finite-element  method. 

A  number  of  numerical  methods  have  been  used  for  solving  heat  transfer  problems  with  phase 
change.  The  essential  problem  is  to  secure  a  solution  for  the  moving  boundary.  The  solution 
methods  are  usually  divided  into  two  basic  types.  The  first  type  is  represented  by  using  tradition¬ 
al  fixed  mesh  finite  elements  or  finite  differences  where  the  phase  front  progresses  through  the 
stationary  mesh  and  is  interpolated  with  temperature  distributions.  Latent  heat  is  handled  with 
the  methods  of  enthalpy  or  apparent  heat  capacity.  In  general,  phase  change  occurs  over  a  very 
narrow  temperature  range  where  enthalpy  is  represented  by  a  steep  line  and  apparent  heat  capac¬ 
ity  has  a  peak  value  (Fig.  1 ;  O'Neill  1 983).  An  artificially  extended  temperature  range  AT  is 
usually  used  to  smooth  the  curves  for  both  parameters.  In  practice,  numerical  results  are  sensi¬ 
tive  to  the  selection  of  AT,  and  any  unreasonable  value  of  AT  will  introduce  physical  distortion. 


O’Neill  (1983)  presented  an  approach  for  the 
case  of  two-dimensional  problems;  with  the 
use  of  linear  triangular  elements,  the  effect 
of  latent  heat  was  included  in  apparent  heat 
capacity,  which  is  theoretically  infinite  at  a 
discrete  phase  change  temperature.  Thus  the 
artificially  extended  temperature  range  AT 
is  not  needed.  This  approach  has  the  advan¬ 
tages  of  simplicity  and  efficiency. 

The  second  type  of  solution  method  in¬ 
volves  the  use  of  a  continuously  deforming 
mesh,  such  that  the  phase  boundary  always 
lies  on  a  particular  numerical  boundary.  Be¬ 
cause  of  the  distinct  boundary  between  the 
two  phases,  every  element  has  its  unique 
physical  properties.  Recently,  Lynch  and 
O’Neill  (1981)  and  O’Neill  and  Lynch  (1981) 
developed  a  method  where  finite  differences 
in  time  are  used  and  the  mesh  motion  effects 
appear  as  a  velocity  term  in  the  governing 
equation.  This  type  of  method  is  considered  Fjgurg  ;  Enthalpy  and  apparent  heat 

to  be  more  flexible  and  accurate  than  the  ity  over  a  temperature  range. 

fixed  mesh  method  and  also  capable  of  solv¬ 
ing  some  special  problems,  e.g.  the  simulation  of  ice  crystal  growth  and  the  process  of  icing  of 
flowing  water  in  a  pipe  (Sullivan  et  al.  1985,  Albert  1984,  Albert  and  O’Neill  1985).  In  this 
paper,  both  fixed  and  deforming  mesh  finite-element  methods  are  used.  In  using  the  deforming 
mesh  finite-element  method,  an  automatic  mesh-generation  technique,  transfinite  mappings 
(Albert  1984)  at  each  time  step  are  adopted. 


BASIC  FINITE-ELEMENT  FORMULAS 


Description  of  problem 

Figure  2  shows  a  cross  section  of  a  cylindrical  shaft,  where  r  is  the  excavation  radius  and  R  the 
radius  of  freeze-pipe-circle.  As  is  shown  in  the  figure,  neither  the  distance  between  adjacent  pipes 
nor  between  the  pipes  and  the  center  is  equal.  Consequently,  the  thickness  of  the  global  frozen 
wall  is  not  uniform.  During  the  active  freezing  stage,  the  global  frozen  wall  will  be  advancing  un¬ 
evenly  in  both  inward  and  outward  directions. 

The  governing  equation  to  be  solved  in  . 


each  phase  is  the  classical  heat  conduction 
equation: 


C  ~  =  V  ■ 
Ld  t  V 


(KVT)  +  Q 


(1) 


with  the  interface  boundary  condition 


d  s 


L  jt  =  {KVT)f-(KVT\ 


(2) 


Figure  2.  Distribution  of  freeze  pipes. 


T  =  temperature 
t  =  time 

C  =  volumetric  heat  capacity 
K  =  heat  conductivity 
L  -  latent  heat 

s  =  position  of  phase  boundary 
Q-  interior  heat 
/  =  frozen 
u  =  unfrozen. 

Finite-element  equation -fixed  mesh 

Using  the  Galerkin  finite-element  method,  we  multiply  eq  1  by  a  weighting  function  Ni  and 
integrate  over  the  whole  domain  A : 


J  NX[C^-V  •  (KVT)-Q\  dA=Q. 


Integration  of  the  second  term  by  parte  yields 


fA  [*, C +  W,  •  KVT] dA-fr  N.,K^dy+JA  QN \ dA 


where  T,  indicates  the  boundary  with  normal  direction  n.  Now  let  the  temperature  be  approxi¬ 
mated  as 


X  ( 5 ] 

M  1  ' 

Here  N j  stands  for  the  interpolating  function,  which  is  chosen  to  be  the  same  as  the  weighting 
function  jV(  . 

Substituting  eq  5  in  eq  4  gives 

jf  C/Vj Nt  dA  +TjfA  VA'j  •  KVNi  dA  -  ^  NyK  dy  +  QNX  dA  =  0.  (6] 

Discarding  the  term  containing  the  integration  of  boundary  heat  flux  and  writing  eq  6  in  matrix 
notation,  it  becomes 


\A  ]  rj  +  [A  ]  T  +  jr  =  0, 


[A]  =  f  CN.NdA 

Ja  ’ 

(A)  =  f  VA'j  •  KVN-dA 
J  A 


M  -L QN id4- 

The  detailed  derivation  of  If  j  is  shown  in  Appendix  A. 


A  key  problem  remaining  to  be  solved  is  the  integration  of  the  integral  that  contains  heat 
capacity  C.  To  take  into  consideration  the  effect  of  latent  heat  L,  O’Neill  (1983)  presented  an 
application  of  a  5  function  by  relating  volumetric  sensible  heat  cs  and  latent  heat  L  as 

C  =  cs+L6  (T-Tp)  (8) 

where  C  is  the  volumetric  apparent  heat  capacity.  It  should  be  noted  that  with  the  use  of  a  5 
function,  phase  change  will  occur  at  a  discrete  temperature  T  . 

Substituting  eq  8  in  eq  6,  the  first  term  gives 

/  c^N^dA  +  f  LNtN}b{T -  Tp)dA.  (9) 

A  A 

The  second  term  in  eq  9  occurs  only  in  the  elements  which  contain  the  Tp  isotherm,  a  detailed 
evaluation  of  which  is  shown  in  Appendix  B.  In  the  case  of  linear  elements,  K  and  cs  may  be 
evaluated  using  the  weighted  area  method. 

Finite-element  equation-deforming  mesh 

The  Galerkin  finite-element  equation  for  deforming  mesh  takes  the  same  form  as  eq  4: 

/|  [v,  r 37  +  VN'  •KVT\dA~fr  NiK  +  [A  QN' dA  =0'  (1 0) 

But  due  to  the  deforming  of  the  mesh  with  time,  the  interpolating  function  N-  is  now  a  function 
of  the  mesh  as  well  as  time  r.  In  a  Cartesian  system,  let 


r*  Z  Tp )  Afj(x,  v,  f ) . 

Substituting  eq  1 1  into  eq  10,  3770;  will  yield  two  terms  and  eq  10  becomes 


r  _  r  3A;  dT.  r 

T  /  AVM  •  VN{dA  +T  CN'-r-l  dA  +-T1  /  CN-.N-.c 
1 J a  >  1  >  Ja  'dr  dt  Ja  '  > 


-  §  N  K  —  dy  +  f  QN-.dA  =  0. 
r,  1  <•»  J  a  ' 


Lynch  (1982)  indicated  that  the  value  of  dN-/dt  is  expressed  in  terms  of  Kby 


9M 

_J=_  V.VN, 
dt  1 


where  V -  (dxjdt)^^  Xj  denotes  the  coordinates  of  nodes  with  respect  to  a  fixed  reference  frame, 
i//j  is  the  inte  polating  function  and  V  is  the  mesh  velocity.  For  the  special  case  of  linear  triangular 
elements,  =  N-v 

Substituting  eq  13  in  eq  12  gives 


T,  f  KVN  •  VN  dA  -  T.  f  CN-  V  •  VN-dA  +~  f  CNN  dA 
1 JA  j  JJA  1  ja  J 

-*rtNiK^dy+fAWidA  =  0-  <14) 

[v4 1 1 7- j  +  [«t]|r|  +  |Fj  =0.  (15) 

Equations  1 5  and  7  are  in  the  same  form,  but  the  difference  between  them  is  in  the  representation 
of  [AT] .  [A'] ,  in  this  case,  is 

[A]  =  |  KVN-VNMA-  CNV.VNMA.  (16) 

JA  '  JA  1 

The  solution  of  the  second  term  in  eq  16  is  shown  in  Appendix  C. 

Equations  14  and  15  are  the  basic  finite-element  equations.  To  specify  the  motion  of  the  phase 
front,  two  methods  are  used  in  this  investigation. 

Method  1 

The  following  derivation  appears  as  it  was  presented  in  Albert  (1984).  Equation  2  is  the  base 
for  calculating  the  motion  of  phase  boundary  nodes.  Since  eq  2  cannot  be  exactly  satisfied  at  a 
discretized  boundary.  Lynch  (1982)  proposed  the  use  of  a  weaker  integral  form  of  eq  2: 

/r2  /.(f]  A'  •  n dy  =/rj  (kVt){  •  n dy  •  n dy 

= /  (a'Vt),  •  nZN^y  -f  (kVT)u  ■  nW^y  (17) 

f  2  1^2 

where  XA,J  =  1 ,  /  refers  to  nodes  on  the  phase  boundary,  n  is  the  unit  vector  directed  away  from 
the  frozen  zone,  and  the  integration  is  over  the  phase  boundary  T2 .  If  we  now  consider  each 
boundary  node /,  eq  17  may  be  written 


'-(jX- L  -vr  ■/.  |(‘  S ),  -  (*  ft  )„ I  W 


where  BT/dn  is  the  temperature  gradient  normal  to  the  local  boundary  and  ( ds/dt  )j  is  the  velocity  of 
node  /  with  the  unit  vector  nij.  By  using  the  characteristics  of  linear  triangular  elements  and  repre¬ 
senting  two  adjacent  sides  of  node  /  with  length  E,  and  E2 ,  we  have 

I  nN^dy  =  'A  Cj  n2  +  Vi  C2  n2  (19) 

■'r2 


,,  __m  __ ,  K4D,-(4"i. 

1  \dt ) ■  /.  Cj  nij  •  n2  +  S2  mj  •  n2 


The  magnitude  of  I  j  may  be  evaluated  from  eq  20  if  mj  is  specified.  To  circumvent  this  difficulty, 
nij  is  assumed  to  be  a  weighted  average  of  n,  and  n2 : 


mj  Sin,  +  S2n, 

‘i^i '  "TOT 


It  should  be  noted  that  the  direction  of  mj  can  be  chosen  at  one’s  convenience;  consequently  the 
magnitude  of  Vi  will  be  in  the  chosen  direction  (App.  D  demonstrates  this  procedure). 


Method  2 

In  a  later  paper,  Lynch  (1983)  pointed  out  the  deficiency  of  method  1  and  reported  a  complete 
heat  conservation  method.  The  refined  basic  Galerkin  equation  is 


I  C  ff  NidA  +f  KVT  •  V/VfdA  +  J  QNidA  ~  F, 


where  F{  includes  two  components,  Ffl  and  Ff*,  which  are  the  integrals  of  heat  flux  over  the  ele¬ 
ment  and  phase  boundaries  respectively,  i.e. 


Ft  =F[>  +  Ff>  =  KVT  •  n N-tdy  +  J  |(aV  rj{  -  (kVT^J-  n N-tdy 


The  node  motion  on  the  phase  boundary  was  evaluated  through  the  combination  of  eq  18, 22,  and 
23  to  form  a  complete  heat  conservation  equation: 


J  C^NidA  +  KVT’  VNdA  +  /  QN,dA 
A  <>'  '  JA  '  JA 

’*r,  *©  WLV>i 


Lynch  (1983)  showed  that  the  new  phase  boundary  condition  yields  second-order  accuracy  and  is 
easy  to  implement  without  extra  computational  difficulties  (App.  E  shows  the  computational  pro¬ 
cedure). 


TRANSFINITE  MAPPING  TECHNIQUE 


In  both  of  the  deforming  mesh  methods  used  here,  the  interior  mesh  motion  is  governed  by  a 
mesh  generation  technique.  Albert  (1984)  reviewed  such  techniques  and  demonstrated  the  use¬ 
fulness  of  transfinite  mappings  (Haber  et  al.  1981)  in  conjunction  with  moving  meshes.  The  inter¬ 
ior  mesh  motion  in  this  work  is  based  on  the  method  presented  by  Albert  (1984). 

In  this  method,  we  introduce  a  concept  of  a  lofting  projector  P  which  maps  a  true  surface  to  an 
approximate  surface  with  a  linear  interpolatory  constraint.  In  Figure  3,  i//j ,  \p2  and  ,  £2  are 
four  curvilinear  boundaries  of  a  plane  region,  and  we  have  a  bilinear  projector: 


P(u,v)  =  (1-vWi  («)+  v\p2  (u)  +  (l  -u)?j  (!’)  +  m£2  (v) 


-  (1  -II) (1  -  v)F(0,0)  -  (1  -u)vF( 0,1) 


-  uvF  (1,1)  -i/(l  -v)F(l,0) 


0  <  u<  1  ,  0  <  r  <  l 


-.  \  .%  .v 


where  u  and  v  are  normalized  coordinates  that  change  linearly  along  ,  'Pi  and  £1 ,  £3,  respective¬ 
ly.  For  a  three-sided  region  containing  curvilinear  boundaries  \p,  £  and  v,  the  trilinear  projector  is 


P(u,v,W)  =  *  (v)  +  (~j  1}  (1-v)  +  (j~)  V  (w) 

+  fe)  *  (1-W)  +  {T^)*  (U)  +  fc)  *  (1-“) 

-  w  (0)  -  u£  (0)  =  vtj  (0) 

0<u<l,0<v<l,0<w<l  and u  +  v  +  w  =  1  .  (26) 


As  before,  u,  v  and  w  are  normalized  coordinates  that  change  linearly  along  \p,  £  and  17. 

Applying  eq  25  and  26,  the  boundaries  of  plane 
region  may  be  divided  into  any  number  of  discrete 
points  as  well  as  unevenly  located.  The  trilinear  pro¬ 
jector  results  in  triangular  elements  directly,  while 
the  bilinear  projector  forms  unit  squares  that  can  be 
transformed  to  triangular  elements  with  diagonals. 

In  principle,  eq  25  and  26  may  be  applied  to  any 
plane  region  of  complex  shape.  This  is  done  by  di¬ 
viding  the  whole  region  into  a  number  of  sub-regions. 

On  occasion,  it  is  more  convenient  to  use  bilinear  and 
trilinear  projectors  to  deal  with  three-  and  four-sided 
regions  respectively.  Albert  (1984)  discusses  the  use 
of  these  projectors  and  provides  guidance  on  situa¬ 
tions  where  they  may  fall. 

At  the  end  of  each  time  step,  a  new  mesh  is  gener¬ 
ated  with  reference  to  the  new  phase  front  and  other 
boundaries,  and  which  forms  the  basis  of  the  compu¬ 
tation  for  the  next  time  step. 


COMPUTATIONS  AND  CONCLUSIONS 

In  this  investigation,  only  a  quarter  of  a  circular 
region  with  ten  freeze-pipes  is  modeled.  The  physi¬ 
cal  parameters  used  here  are  obtained  from  the  Figure  3.  Mesh  generation  with  trans- 

Panji-3  Dong  Feng  Shaft  of  the  Liang-Huai  Coal  finite  mappings. 

Mine  in  eastern  China  (Table  1). 

The  two  straight  boundaries  of  the  computational  region  are  specified  to  be  zero-flux  boundary 
conditions.  In  deforming  mesh  approach,  the  outer  boundary  nodes  are  kept  at  a  constant  temper 
ature  and  moved  outward  at  each  time  step.  The  extent  of  this  movement  is  always  twice  that  of 
the  phase  boundary. 

Three  programs  were  used:  FEFIX  performs  the  fixed  mesh  finite-element  calculations  as  pre¬ 
sented  by  O’Neill,  MOVGR  uses  the  moving  mesh  with  transfinite  mappings  as  presented  by  Al¬ 
bert  (1984),  and  MOVHE  incorporates  the  improvement  over  Albert’s  method  suggested  by 
Lynch.  The  flow  charts  of  the  programs  are  shown  in  Appendix  F.  For  FEFIX,  MOVHE  and 
MOVGR,  there  are  42,  20  and  13  time  steps  adapted  to  simulate  a  period  of  330  days.  The  pro¬ 
gram  MOVGR  iterates  once  on  the  location  of  phase  front  for  each  time  step.  In  all  the  three 
programs,  the  parameter  for  implicitness  in  time  6  is  equal  to  1 .0;  i.e.  it  is  fully  implicit. 


<1.1 ) 


,<»> 


(O.I.O) 


(0,0,1) 


(1.0.0) 


Table  1.  Parameters  used  in  calculations. 


Initial  temperature, 

Boundary  temperature, 

Phase  change  temperature, 

1-feat  conductivity, 

Heat  capacity, 

I 

Latent  heat, 

Dry  density, 

Relative  water  oontent, 

Relative  unfrozen  water  oontent,  I 
Heat  flux  on  pipe  wall. 

Diameter  of  freeze-pipe, 

Statistical  coefficients, 


15°C 

15°C 

0°C 

1.00  kcal/m  •  hr  •  °C 
1.50  kcal/m  •  hr  •  °C 
710  kcal/m1  •  “C 
560  kcal/m*  -°C 
80./XB'-iyu) 

1360  kg/m’ 

26.5% 

5% 

QT-a  f*.  Qr=  250  n d 

0.127  m 

4.08 

-0.28 


Figure  4  shows  the  structure  of  fixed  mesh 
and  Figures  5, 6  and  7  illustrate  the  locations 
of  phase  fronts  at  90, 210  and  330  days  re¬ 
spectively.  Similarly,  for  the  deforming  mesh 
used  with  method  1  and  method  2,  the  phase 
boundaries  for  identical  periods  of  90,  210 
and  330  days  are  shown  in  Figures  8, 9  and 
10  and  Figures  11,12  and  13,  respectively. 
Figure  14  shows  the  comparison  of  thickness 
of  inner  and  outer  frozen  wall  as  a  function 
of  duration  of  operation  among  the  computa¬ 
tional  methods.  Figure  15  shows  the  growth 
of  total  frozen  wall  thicknesses  as  a  function 
of  time.  The  numerical  results  are  also  shown 
in  Table  2.  All  the  results  in  Figures  14  and 
15  and  in  Table  2  are  average  values  of  frozen 
thickness  at  each  time  step. 

In  practice,  the  Fixed  mesh  method  is  found 


Figure  4.  Fixed  mesh  and  freeze-pipe  distribu¬ 
tion  (FEFIX). 


to  be  easier  to  program  and  takes  much  less  CPU  time  than  the  deforming  mesh  method  during  the 
run.  It  is  unfortunate  that,  because  of  the  absence  of  experimental  data  and  analytical  solution, 
the  computed  results  can  only  be  compared  with  those  reported  by  Ding  et  al.  (1982)  using  the 
finite  difference  method  and  under  the  same  physical  conditions.  The  results  of  Ding  have  been 
compared  with  the  field  results  and  have  been  proved  to  be  in  good  agreement. 


Time  (90  days) 


Table  2.  Computed  thicknesses  of  frozen  wall. 

ysj  Time  (I SO  days)  Time  (210  days)  Time  (270  days)  Time  (S00  days) 


Thickness  (m) 

Thickness  (m) 

Thickness  (m  ) 

Thickness  (m ) 

Thickness  (m) 

Inner  Outer  Total 

Inner  Outer  Total 

Inner  Outer  Total 

Inner  Outer  Total 

Inner  Outer  Total 

FEFIX 

1.26 

1.12 

2.38 

1.98 

1.67 

3.65 

2.71 

2.13 

4.84 

3.37 

2.52 

5.89 

3.80 

2.74 

6.54 

MOVGR 

1.36 

1.20 

2.56 

2.05 

1.78 

3.83 

2.69 

2.28 

4.97 

3.31 

2.74 

6.05 

3.62 

2.95 

6.57 

MOVHE 

1.34 

1.19 

2.53 

2.04 

1.75 

3.79 

2.70 

2.23 

4.93 

3.34 

2.66 

6.00 

3.65 

2.86 

6.51 

FINITE 

DIFFERENCE 

1.50 

1.36 

2.86 

2.10 

1.83 

3.93 

2.70 

2.28 

4.98 

3.30 

2.65 

5.95 

3.60 

2.82 

6.42 

8 


sum 


Im 


Figure  8.  Deforming  mesh  at  90  days 
( MOVGR ). 


Figure  9.  Deforming  mesh  at  210  davs 
(MOVGR). 
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Figure  10.  Deforming  mesh  at  330  days 
(MOVGR). 


Figure  11.  Deforming  mesh  at  90  days  (MOVHE). 
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Figure  12.  Deforming  mesh  at  210  days 
(MOVHE). 


Figure  13.  Deforming  mesh  at  330  days 
(MOVHE). 


360 

1 - 1 - 

1 - 1 - 1 - 1 - 1 - 1 

— 

- 

0* 

0  FEFIX  Ofe 

— 

300 

_ 

A  MOVGR 

_ 

o 

•  MOVHE  0 

— 

A 

A 

- 

240 

— 

• 

• 

— 

- 

% 

0* 

- 

180 

— 

• 

• 

— 

- 

- 

120 

0  o 

— 

— 

z*  ^ 

— 

60 

— 

£q 

— 

INNER 

•  •  OUTER 

— 

<£ 

— 

_ 1 

1  1  1  1  1 

-  1 

f 

Jk 

1 

IM 

-2-10123 
Frozen  Thickness  (m) 

4 

Figure  14.  Changes  in  inner  and  outer  thickness. 


Figure  15.  Changes  in  total  frozen  thickness. 
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APPENDIX  A:  POINT  HEAT  SOURCES. 


In  eq  7  and  1 5 , 

\f\-Ia  QNjJA  (Al) 

For  elements  containing  freeze-pipe,  the  heat  source  Q  is  treated  as  a  point  source  and  by  introduc 
ing5  functions  (Segerlind  1984),  eq  1  becomes 

f  QNfi  (x  -x0)6  O’  - }’0 ) dxtJy  ( A2 ) 

J  A 


where  (x0,  v0)  is  the  coordinate  of  freeze-pipe.  Thus  we  have 


*i(Wo) 
^  Ar  (-V0 ,  v0 ) 


(A3) 


where  i,  j,  k  denote  the  three  nodes  of  the  element.  By  using  eq  A3,  the  effect  of  Q  may  be  auto¬ 
matically  distributed  to  neighboring  nodes,  depending  on  the  relative  location  from  the  point 


APPENDIX  B:  EVALUATION  OF  THE  INTEGRAL  INCLUDING  LATENT  HEAT  (FIXED 
MESH). 


The  integral  including  latent  heat  is 

f  LNppW-TJdxdy.  (Bl] 

To  integrate  this,  one  of  the  spatial  variables  must  be  expressed  in  terms  of  T.  In  Figure  Bl ,  the 
isotherm  crosses  the  element  and  intersects  with  two  of  the  sides  at  points  1  and  2.  Let  the  new 
coordinate  system  be  ( x',y ),  where  y  extends  along  the  isotherm,  and  x  is  perpendicular  to  the 
isotherm.  Equation  1  becomes 


fim 


6  ( T-Tp)dx'dy ' . 


where  V;  and  Nj  are  the  functions  of  x  and.v’  Thus  all  isotherms  in  the  element  are  parallel  to 
y  ,  T  is  only  a  function  of  x  and  dx'  is  ( dx’IdT)  •  dT.  Therefore  eq  B2  can  be  written  as 


where  dx'jdT,  the  reciprocal  of  temperature  gradient  along  x',  is  a  constant.  Applying  the  charac 
ter  of  5  function,  eq  B3  becomes 

ry 2 


f  Vj/V  I 

Jy\  '  Mr-: 


Equation  B4  is  a  line  integral  along  the  isotherm.  Using  Simpson’s  rule,  eq  B4  becomes 

Ldf’  !  +^«a^i.) 


-fo=(biTi+bjTj+bkTii)/2A 


3 7m<CtTi  +  clTi+ckTkV2A. 


fi  represents  the  length  of  isotherm,  which  can  be  evaluated  geometrically.  The  values  of  A'. , ,  Ni2 , 
jVj,  and  Afj  can  easily  be  specified  because  of  the  properties  of  interpolation  functions,  and  in 
eq  B6,  i,  j  and  k  denote  element  nodes. 


Figure  Bl.  Evaluation  of  the  integral  including 
latent  heat. 
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APPENDIX  C:  SPECIFICATION  OF  (X]  -DEFORMING  MESH. 

There  are  two  terms  in  [X] .  One  of  the  terms  is 

/  CN-.V'V NdA  (Cl 

JA  ’ 

where  C  is  constant  for  any  element  at  each  time  step,  and  Vm  VAj  may  be  written  as  Vs  ctNJbx^ 

=  j_  [b.t  bj  bkl  (c2 

dxe  2A  c-,  c}  rk 

where  6  denotes  x  and_y,  and  i,j  and  k  denote  three  nodes.  Ke  can  be  written  as 


V  , 

V  . 

XI 

y* 

V*> 

*yi 

^yk 

where  Kxi,  Vyi  and  so  on  can  be  calculated  from  node  locations  for  the  two  meshes  representing 
the  region  at  the  beginning  and  end  of  the  time  step.  Using  eq  C2  and  C3  in  eq  Cl ,  we  have 

CV.VN.  (  N-N.dA  .  (C4 


V 


APPENDIX  D:  SPECIFYING  dT/dn  AND  THE  DIRECTION  OF  rrij  FOR  METHOD  1 . 

When  the  radial  direction  of  nij  is  specified,  the  solution  of  •  nt  and  £2nij  *  n2  may  be 
implemented  geometrically.  In  Figure  D1 , 

mj*ni  =  |mj|  •  Injl  •  co sfi  =  di/d2 

Ejm^n,  =  x^y^/^-xf  +0'j-7j)2  •  (E 

£2nyn2  can  be  evaluated  in  the  same  manner. 


\  dJ  / 

d,\  / 

V/ 


Figure  Dl.  Specification  of  the  motion  of 
node  j  on  phase  boundary. 


dT/dn  is  the  temperature  gradient  normal  to  phase  boundary.  In  eq  2  (in  text),  dT/dn  contains 
four  terms  referring  to  the  adjacent  two  sides  and  to  the  frozen  and  unfrozen  zones  respectively. 
Each  may  be  specified  with  the  following  equation: 


dT  fojsr  /drv2 
dn  \l\dx)  \3y/  1 


APPENDIX  E:  PROCEDURES  OF  METHOD  2. 


The  complete  heat  conservation  equation  (eq  24)  (in  text)  may  be  written  in  matrix  notation 


I  !  _  2  (Af]  M"*1  -  (A 
I  M  L  C,n,  +  £:n2 


(El 


where 


[Af]  =  [A]  +  0Ar[/t] 

|iv|  =  ([A]  -(1  -6)At[k})  rj"  +  At[(l  - 
The  procedure  for  computing  Is,  I  is  as  follows: 


1 ,  Form  global  equations  as  usual.  .  . 

2,  Save  the  equations  of  phase  boundary  nodes,  i.e.  [Af  ]  and  \N )  . 

3,  Solve  the  nodal  temperature. 

4,  Compute  the  latent  heat  balance  with  the  saved  equations  in  step  2  and  the  new  temperature 
found  in  step  3. 


appendix  f:  explanation  of  programs. 


(A)  Name 

FEFIX  -  Fixed  mesh  finite  element  method  (flow  chart  in  Figure  FI ), 
MOVGR  -  Deforming  mesh  with  method  1  (flow  chart  in  Figure  F2). 
MOVHE  —  Deforming  mesh  with  method  2  (flow  chart  in  Figure  F3). 


Figure  FI.  Flow  chart  for  FEFIX. 


Figure  F2.  Flow  chart  for  MOVGR.  Figure  F3.  Flow  chart  for  MOVHE. 


(B)  FEFIX 

FEFIX  contains  the  main  program  and  subroutines  GLOB,  MAT,  PHASE,  RATIO  and  BAND. 
GLOB  -  Forms  global  matrices  and  inserts  boundary  conditions. 

MAT  —  Specifies  the  element  matrices. 

PHASE  -  Specifies  the  location  of  the  T  isotherm. 

RATIO  —  Specifies  parameters  A  and  C  in  the  elements  containing  phase  change. 

BAND  -  Solves  symmetric,  banded  and  positive  definite  matrix  equations  and  obtains  nodal 
temperatures. 


(C)  MOVGR 

MOVGR  contains  main  program  and  subroutines  MOVE,  MESH,  AREAS,  GLOB,  MAT  and 
BAND. 

MOVE  --  Specifies  the  motions  of  phase  boundary  nodes  with  method  1. 

MESH  -  Generates  a  new  mesh  using  the  new  locations  of  phase  boundaries  and  transfinite 
mappings.  The  velocities  of  nodal  motion  are  calculated  from  the  locations  of  new 
and  old  meshes. 

AREAS  —  Specifies  new  element  areas  at  each  time  step. 

GLOB  -  Forms  global  matrices  and  inserts  boundary  conditions. 

MAT  —  Specifies  the  element  matrices. 

BAND  —  Solves  banded  and  positive  definite  matrix  equations  and  obtains  nodal  temperatures. 

(D)  MOVHE 

MOVHE  contains  main  program  and  subroutines  MOVE,  MESH,  AREAS,  GLOB,  MAT  and 
BAND.  In  the  subroutines,  only  MOVE  and  GLOB  are  different  from  those  in  MOVGR. 

MOVE  —  Specifies  the  motions  of  phase  boundary  nodes  with  method  2,  in  which  the  saved 
nodal  equations  are  used. 

GLOB  -  Forms  global  matrices,  inserts  boundary  conditions  and  saves  the  nodal  equations 
on  phase  boundary. 


A  facsimile  catalog  card  in  Library  of  Congress  MARC 
format  is  reproduced  below. 
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